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Preface 

The  work  presented  in  this  report  was  performed  at  Matis,  Inc.,  1255  Biltmore  Drive, 
Atlanta,  Georgia  30329  and  at  the  Department  of  Mathematics,  University  of  California, 
Los  Angeles,  CA  90095-1555.  This  work  was  sponsored  by  the  Air  Force  Office  of  Scientific 
Research  during  the  period  September  1,  1998  -  August  31,  2000.  The  project  Technical 
Monitor  was  Dr.  Arje  Nachman  from  the  AFOSR.  We  are  indebted  to  Dr.  Nachman  for 
all  the  encouragement  and  support  he  gave  us  during  the  work  on  this  contract. 
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1  Executive  Summary 

Two  directions  of  work  have  been  pursued  under  this  effort.  The  first  direction  is  concerned 
with  further  extensions  of  our  geometric  techniques  for  spatial  and  surface  ray  tracing 
and  applications  of  these  techniques  to  problems  in  high-frequency  electromagnetics.  In 
addition,  high-frequency  physical  optics  methods  were  developed  for  accurate  and  efficient 
calculation  of  scattered  fields  due  to  currents  in  the  penumbra  region.  The  second  direction 
is  concerned  with  development  of  grid  methods  for  capturing  the  motion  of  self-intersecting 
interfaces  and  related  PDEs  connected  with  propagation  of  high-frequency  waves.  In  both 
directions  significant  progress  has  been  achieved.  For  convenience  of  reference  these  two 
directions  are  designated  below  as  Efforts  A  and  B. 

Under  effort  A,  the  theoretical  work  carried  out  by  the  Matis,  Inc.  team  resulted  in 
new  methods  for  solving  global  problems  of  determining  geodesics  in  presence  of  obstacles 
with  highly  complex  geometries.  This  is  a  critical  capability  necessary  for  finding  accurate 
and  reliable  high-frequency  solutions  to  numerous  scattering  problems.  The  theoretical 
work  was  accompanied  with  a  very  extensive  developments  of  computational  methods. 

In  particular,  the  theoretical  results  obtained  under  effort  A  were  implemented  into 
a  system  of  computational  codes  integrated  into  two  fully  operational  computer  software 
systems  for  interactive  analysis  of  airborne  antenna-to-antenna  electromagnetic  interfer¬ 
ence/compatibility  and  interactive  analysis  of  airborne  antenna  radiation  patterns.  Both 
computer  systems  have  truly  revolutionary  capabilities  to  perform  antenna  analysis  on  real¬ 
istic  digital  aircraft  models  represented  by  models  produced  with  commonly  used  computer 
aided  design  (CAD)  methods.  Even  though  the  development  of  these  systems  has  only  re¬ 
cently  been  completed,  a  number  of  inquiries  have  been  received  from  the  Air  Force  and 
commercial  companies  working  on  defense  contracts  with  the  US  Government.  Such  inter¬ 
est  in  the  developed  systems  is  strong  evidence  that  the  results  of  this  work  are  applicable 
to  solution  of  real-life  problems  important  for  the  US  defense. 

In  addition,  under  the  effort  A  the  Matis  team  developed  a  very  efficient  and  stable 
method  for  computing  the  correction  terms  to  the  physical  optics  (PO)  fields  due  to  currents 
at  the  shadow  boundary  of  a  scatterer.  This  method  will  be  applied  to  enhance  the  accuracy 
of  radar  cross  sections  (RCS)  predictions  by  PO  methods. 

Under  effort  B,  the  team  at  UCLA  developed  a  revolutionary  fixed  grid  method  for 
capturing  the  motion  of  self-intersecting  interfaces  and  related  PDEs.  Moving  interfaces 
that  self-intersect  arise  naturally  in  the  geometric  optics  model  of  wavefront  motion.  Ray 
tracing  techniques  can  be  used  to  compute  these  motions,  but  they  lose  resolution  as 
rays  diverge.  This  difficulty  was  overcome  by  developing  a  new  numerical  method  that 
maintains  uniform  spatial  resolution  of  the  front  at  all  times.  The  approach  is  a  fixed 
grid,  interface  capturing  formulation  based  on  the  Dynamic  Surface  Extension  method 
of  Steinhoff  and  Fan.  The  new  methods  can  treat  arbitrarily  complicated  self-intersecting 
fronts,  as  well  as  refraction,  reflection  and  focusing.  We  also  further  extended  this  approach 
to  curvature  dependent  front  motions,  and  the  motion  of  filaments.  The  new  methods  have 
been  validated  with  numerical  experiments. 

In  summary,  both  directions  of  work  were  complementary  to  each  other  and  produced 
significant  theoretical  results  and  numerical  codes  with  capabilities  for  solving  important 
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and  difficult  high-frequency  wave  propagation  problems. 
We  describe  now  the  results  in  more  detail. 


2  Effort  A 

2.1  Geometric  Techniques  in  High-Frequency  Computational  Elec¬ 
tromagnetics 

An  important  practical  problem  in  communications  is  to  numerically  predict  radiation 
directivity  properties  of  antennas  mounted  on  geometrically  complex  platforms  acting  as 
scatterers.  High-frequency  methods  of  electromagnetics  have  been  used  extensively  for 
dealing  with  this  problem  in  cases  when  the  platform  can  be  represented  as  a  collection 
of  a  small  number  of  simple  shapes  such  as  cylinders,  ellipsoids,  pieces  of  planes,  cones, 
etc.  When  such  a  representation  of  the  platform  is  sufficiently  accurate  the  high-frequency 
methods  provide  an  accurate  estimate  of  the  electric  field  produced  by  an  antenna  in  a 
given  direction.  However,  comparisons  of  computed  results  with  experimental  data  showed 
that  large  errors  may  result  in  from  inaccurate  representation  of  the  platform  geometry. 
Consequently,  there  is  a  need  for  algorithms  capable  of  performing  high-frequency  analyses 
in  geometrically  complex  environments  that  can  not  be  reduced  to  small  collections  of 
simple  shapes. 

High-frequency  ray  methods  are  based  on  geometrical  optics  and  its  extensions  and,  in 
order  to  apply  these  methods,  one  needs  first  to  construct  the  propagation  paths  from  the  ' 
antenna  to  a  given  observation  point  at  infinity  specified  as  a  point  on  a  unit  far-sphere*. 
Such  observation  point  at  infinity  is  usually  referred  to  as  the  “far-field”  direction.  When 
the  scatterer  does  not  obstruct  the  wave  propagation  the  path  is  a  linear  ray.  Otherwise, 
the  scatterer  acts  as  an  obstacle  and  the  propagation  path  may  have  a  nonempty  contact 
set  with  the  surface  of  the  scatterer.  In  the  present  setting  the  scatterer  is  assumed  to  be 
perfectly  conducting,  and  no  penetration  of  the  paths  through  the  scatterer  is  allowed. 

One  of  the  main  principles  of  the  classical  geometrical  optics  is  that  fields  propagate 
along  paths  that  are  critical  points  of  the  Fermat  functional.  The  Euler  equations  for  these 
critical  points  lead  to  the  eikonal  equation.  In  classical  setting  the  scatterer  is  smooth  and 
usually  only  the  paths  that  propagate  in  free  space  and  reflect  or  refract  are  considered. 

In  contrast,  in  the  geometric  theory  of  diffraction  and  its  extensions  (especially,  the  uni¬ 
form  theory  of  diffraction)  the  notion  of  a  propagation  path  is  generalized  to  include  paths 
not  only  with  spatial  segments  and  reflections,  but  also  those  that  diffract  at  geometric 
singularities  such  as  corners  and  edges,  may  attach  to  the  surface  of  the  scatterer,  “creep” 
along  it,  then  detach,  diffract,  and  so  on. 

If  a  propagation  path  from  the  source  to  an  observation  point  is  known,  the  field  at 
the  observation  point  due  to  that  path  can  be  obtained  by  “gluing”  together  the  local 
scattering  information  along  the  path.  In  principle,  this  scattering  information  can  be 
computed  by  high-frequency  techniques,  though  application  of  these  techniques  is  by  no 

’Typically,  an  engineer  may  specify  a  set  of  sample  directions  (“pattern  cut”)  in  which  the  electric  field 
produced  by  an  antenna  is  required. 
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means  straight-forward.  Finally,  the  total  field  at  the  observation  point  is  obtained  as  a 
vector  sum  of  contributions  by  the  fields  propagating  along  individual  paths  from  the  source 
to  that  observation  point. 

The  geometrical  optics  theory  initiated  by  J.  Keller  and  developed  continuously  since 
1950  focused  primarily  on  local  asymptotics  of  diffracted  fields,  without  providing  tools 
for  identifying  and  computing(!)  propagation  paths  in  global,  that  is,  on  entire  platforms. 
In  present  work,  a  theory  of  globally  defined  critical  paths  of  the  Fermat  functional  was 
developed  in  the  case  of  polyhedral  surfaces.  This  case  is  important  because  in  practice 
scatterers  are  often  represented  in  this  form  and,  furthermore,  for  such  scatterers  we  suc¬ 
ceeded  in  constructing  efficient  computational  methods  for  determining  these  paths. 

The  corresponding  variational  theory  associated  with  the  Fermat  functional  had  to  be 
re-examined  and  extended  to  include  non-classical  cases  required  by  diffraction  theory.  In 
certain  cases,  in  order  to  capture  paths  with  specific  properties  the  variations  may  have 
to  be  restricted.  For  example,  in  order  to  obtain  paths  that  contain  a  reflection  point  on 
the  scatterer  one  should  require  that  admissible  paths  have  a  non-empty  contact  set  on 
the  scatterer.  Similarly,  in  order  to  obtains  paths  diffracting  at  corners  or,  in  some  special 
cases  at  edges,  one  needs  to  further  restrict  classes  of  admissible  variations.  The  definitions 
of  stationary  paths  in  such  cases  must  be  modified  accordingly. 

In  applications  related  to  directivity  properties  of  antennas  the  diffracted  field  along 
creeping  segments  of  a  propagation  path  within  the  shadow  region  decays  exponentially. 
Therefore,  it  is  important  to  construct  only  the  “dominant”  paths  and  avoid  searching  for 
paths  with  contributions  to  the  electric  field  below  some  user-specified  level. 

Under  these  effort,  algorithms  for  building  the  dominant  paths  in  the  discrete  radiation 
problem  have  been  developed  and  implemented.  The  algorithm  starts  with  building  a 
possibly  large  set  of  initial  paths  that  are  quasi-geodesics  in  the  metric  of  the  scatterer. 
The  initial  paths  are  terminated  at  the  shadow  boundary  (relative  to  the  current  far-field 
direction).  For  a  given  far-field  direction  the  algorithm  uses  the  local  optimality  condition 
at  the  shadow  boundary  to  filter  out  some  of  the  ’‘nitial  paths.  In  “generic”  cases  this 
procedure  filters  out  most  of  the  initial  paths.  The  initial  paths  that  remain  are  optimized. 
The  optimization  of  each  of  the  paths  uses  a  a  variant  of  unfolding  strategy  combined  with 
lifting  of  paths  off  the  surface  of  the  scatterer  whenever  possible. 

Extensive  numerical  validations  of  the  codes  were  performed  and  the  results  were  com¬ 
pared  with  available  measured  data.  Very  good  agreement  between  the  numerical  results 
and  the  measured  data  has  been  established. 


2.2  Methods  for  Calculating  Scattering  Due  to  Currents  at  Shadow 
Boundaries 

In  typical  physical  optics  (PO)  approximations  used  for  computing  the  scattered  far  field 
the  surface  of  the  scatterer  is  divided  into  the  lit  and  shadow  regions.  The  PO  current 
induced  on  the  scatterer  is  calculated  as 


JPO  = 


2 n  x  Hl  in  the  lit  region 
0  in  the  shadow  region, 
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where  n  is  the  normal  field  on  the  scatterer  and  H1  is  the  incident  magnetic  field.  The 
fact  that  the  JPO  =  0  in  the  shadow  region  impacts  negatively  the  accuracy  of  calculated 
radiated  fields,  since  the  contribution  from  the  current  in  the  shadow  region  is  not  taking 
into  account  the  creeping  waves  propagating  into  that  region.  Neglecting  this  contribution 
may  lead  to  serious  errors  especially  in  cases  where  bistatic  radar  cross-sections  (RCS)  are 
required. 

The  ILDC  (incremental  length  diffraction  coefficients)  method  has  emerged  during  the 
last  two  decades  as  a  technique  that  in  certain  cases  permits  calculation  of  of  a  correction 
to  PO  field  accounting  for  edge  diffraction  for  arbitrary  incidence  and  scattering  angles. 
However,  the  edge  diffraction  problem  has  a  very  specific  structure  that  makes  it  particu¬ 
larly  suitable  for  treating  by  the  ILDC  method.  Nevertheless,  it  seems  natural  to  examine 
the  applicability  of  the  ILDC  method  to  the  case  when  the  shadow  boundary  on  a  surface 
is  treated  as  a  diffraction  line.  This  approach  has  been  analyzed  by  A.  Yaghjian,  R.  Shore, 
and  M.  Woodworth  for  the  case  of  the  sphere  and,  more  recently,  by  A.  Yaghjian  and  R. 
Shore  for  the  case  of  an  ellipsoid  of  revolution. 

The  objective  of  our  work  was  to  extend  the  ILDC  method  to  the  case  of  arbitrary 
convex  surfaces  represented  by  CAD  (computer  aided  design)  -  generated  models.  Several 
serious  difficulties  have  to  be  overcome  when  dealing  with  general  surfaces.  First  of  all,  in 
the  case  of  an  edge  (or  an  ellipsoid),  the  edge  (shadow  boundary  in  case  of  the  ellipsoid)  is 
well  defined,  geometrically  simple,  and  can  be  described  analytically.  On  general  surfaces 
there  is  no  simple  description  of  the  shadow  line  and  such  a  description  has  to  be  built 
by  numerical  means.  Furthermore,  the  geometric  information  regarding  the  surface  in  the 
shadow  zone  also  has  to  be  obtained  numerically.  Secondly,  in  contrast  with  the  edge 
diffraction  problem,  in  the  case  of  the  diffraction  by  shadow  boundary  the  starting  points 
of  the  diffracted  rays  occupy  the  whole  area  of  the  shadow  region  reached  by  the  creeping 
waves.  Thirdly,  in  the  case  of  edge  diffraction  there  are  procedures  that  allow  to  reduce 
the  computation  of  the  radiation  integral  to  a  line  integral.  Such  procedure,  in  general,  is 
not  available  for  the  shadow  boundary  problem  and  needs  to  be  developed. 

In  our  approach,  developed  jointly  with  Dr.  Arie  Michaeli  from  Israel  ve  have  been 
successful  in  recasting  the  ILDC  method  in  a  form  independent  of  a  specific  surface  (such 
as  sphere  or  ellipsoid)  which  makes  the  method  applicable  to  arbitrary  convex  surfaces. 
Further,  suitable  procedures  for  determining  the  ILDC’s  from  far  fields  scattered  by  canon¬ 
ical  problems  have  been  developed.  We  have  now  developed  an  algorithm  that  can  run  on 
fairly  general  convex  surfaces  and  we  are  testing  it  in  various  special  cases.  The  first  results 
look  very  good  and  are  encouraging.  Once  this  stage  is  complete  and  stable  performance 
of  the  algorithms  is  achieved,  we  can  start  working  on  building  codes  which  can  work  on 
models  represented  by  faceted  and/or  NURBS  (Non-Uniform  Rational  B-Splines)  surfaces. 
The  machinery  that  we  built  for  doing  geometric  computations  during  Phase  I  and  our 
other  work  on  scattering  problems  is  a  key  to  this  development.  An  important  application 
for  this  work  is  in  the  area  of  RCS  predictions.  In  particular,  the  current  version  of  the 
widely  used  RCS  code  X-patch  does  not  have  the  capability  to  account  for  the  fields  due 
to  creeping  waves  in  the  shadow  region.  We  plan  to  explore  possible  options  for  adding  on 
such  a  capability  to  the  X-patch  code  with  the  current  developers  of  X-patch. 
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3  Effort  B 

3.1  A  Fixed  Grid  Method  for  Capturing  The  Motion  of  Self- 
Intersecting  Interfaces  and  Related  PDEs 

We  have  developed  a  revolutionary  fixed  grid  method  for  capturing  the  motion  of  self- 
intersecting  interfaces  and  related  PDEs.  We  outline  the  method  and  results  here. 

Moving  interfaces  that  self-intersect  arise  naturally  in  the  geometric  optics  model  of 
wavefront  motion.  Ray  tracing  techniques  can  be  used  to  compute  these  motions,  but  they 
lose  resolution  as  rays  diverge.  We  have  developed  a  new  numerical  method  that  maintains 
uniform  spatial  resolution  of  the  front  at  all  times.  Our  approach  is  a  fixed  grid,  interface 
capturing  formulation  based  on  the  Dynamic  Surface  Extension  method  of  Steinhoff  and 
Fan,  “Eulerian  computation  of  evolving  surface  curves  and  discontinuous  fields”,  Technical 
Report  UTSI,  Tullahoma,  TN.  The  new  methods  can  treat  arbitrarily  complicated  self 
intersecting  fronts,  as  well  as  refraction,  reflection  and  focusing.  We  also  further  extended 
this  approach  to  curvature  dependent  front  motions,  and  the  motion  of  filaments.  We 
validated  the  new  methods  with  numerical  experiments. 

In  the  limit  of  short  wavelengths,  it  is  well  known  that  a  wavefront  moving  through 
a  medium  can  be  described  as  a  moving  surface  with  a  normal  velocity  that  depends  on 
position, 

v  —  c(x)h 

where  n  is  the  local  normal  to  the  front  and  c(x)  is  the  local  wave  speed.  Notable  examples 
include  the  short  wavelength  approximation  of  seismic  and  electromagnetic  pulses,  as  well 
as  the  familiar  example  of  ripples  moving  on  the  surface  of  a  pond.  An  important  feature 
of  this  idealized  wavefront  motion  is  that  intersecting  wavefronts  pass  through  each  other, 
and  also  that  they  reflect  and  refract  off  boundaries. 

Many  interesting  numerical  methods  have  been  developed  to  compute  these  complex 
motions.  The  most  detailed  approach  is  to  discretize  the  governing  wave  equations  directly. 
Unfortunately,  this  approach  is  often  impractical  because  it  requires  that  the  discretization 
resolve  the  short  wavelengths,  which  may  be  thousands  of  times  smaller  than  the  length 
scale  of  interest. 

At  the  other  extreme,  ray  tracing  can  be  used  to  evolve  waterfronts  according  to  geo¬ 
metrical  optics,  here,  the  front  is  represented  using  a  number  of  markers  which  are  moved 
independently.  This  approach  has  the  advantage  of  simplicity,  but  the  markers  tend  to 
diverge  which  leads  to  loss  of  resolution  and  aliasing  of  the  front. 

To  maintain  a  uniform  resolution  of  the  interface,  it  is  natural  to  consider  a  fixed  grid, 
interface  capturing  formulation  such  as  the  Level  Set  method.  Here,  the  wavefront  is 
represented  as  the  zero  contour  of  a  smooth  function  </>,  which  in  turn  evolves  according  to 
the  level  set  equation 

<f>t  +  c(x)\V  <f>\  =  0. 

This  can  be  solved  accurately  and  efficiently  using  numerical  Partial  Differential  Equation 
(PDE)  techniques.  Unfortunately,  the  basic  level  set  method  is  inappropriate  for  treating 
evolving  wavefronts  because  the  solutions  to  this  PDE  will  have  fronts  merge  upon  colliding, 
rather  than  pass  through  one  another. 
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To  obtain  a  fixed  grid  method  appropriate  for  capturing  wavefront  self-intersection 
Steinhoff  and  Fan  proposed  Dynamic  Surface  Extension  (DSE)  methods.  These  schemes 
start  from  some  spatially  distributed  representation  of  the  interface  (similar  to,  but  more 
general  than,  the  level  set  <j>  representation),  and  the  motion  is  achieved  by  alternating 
between  two  steps:  a  simple  short  time  evolution  comparable  to  ray  tracing,  and  an  ex¬ 
tension  step  that  updates  the  distributed  representation  to  reflect  the  new  front  location. 
DSE  methods  automatically  give  a  uniform  resolution  of  expanding  fronts  by  using  a  fixed 
spatial  grid,  and  the  fronts  automatically  pass  through  one  another  rather  than  merging. 

The  original  DSE  methods  are  not  well  suited  to  certain  fundamental  self-intersection 
problems  such  as  the  formation  of  swallowtails.  We  have  generalized  a  basic  DSE  scheme 
(the  Closest  Point  Method)  to  handle  this  fundamental  problem,  as  well  as  all  other  complex 
intersections.  We  further  generalized  our  approach  to  reflecting  and  refracting  wavefronts. 
We  also  discussed  new  extensions  for  propagating  intensity  values,  for  treating  curvature- 
dependent  flows,  and  for  treating  the  motion  of  filaments  (or  more  generally,  objects  of 
co-dimension  >  1). 

The  Closest  Point  Method 

Initialize.  For  each  point  x  €  Rn:  Set  the  initial  tracked  point  TP(x)  equal  to  the  closest 
point  (to  x)  on  the  initial  surface  A0.  set  h  equal  to  the  surface  normal  at  the  tracked  point 
TP(x ),  and  let  c  denote  the  wavefront  speed  at  the  tracked  point. 

Repeat  for  all  steps: 

1.  Evolve  the  tracked  point  TP(x)  according  to  the  local  dynamics  for  a  time  At  : 
TP(x)t  —  cn. 

2.  Extend  the  surface  representation  by  setting  each  tracked  point  TP(x)  equal  to  the 
true  closest  point  on  the  updated  surface  Gamma,  where  T  is  defined  to  be  the  locus 
all  tracked  points,  i.e.,  T  =  {TP(x)\x  e  Rn}.  Replace  each  h(x)  by  the  normal  at 
the  updated  TP(x). 

End. 

Intuitively,  the  manner  in  which  this  method  treats  self-intersection  is  most  easily  under¬ 
stood  by  considering  how  it  treats  two  colliding,  planar  waves.  Initially,  each  nodal  tracked 
point  value  is  set  equal  to  the  closest  point  on  the  nearest  wavefront.  These  tracked  points 
are  updated  during  the  Evolution  Step  according  to  TP(x)new  =  T P{x)or%9mal  -1-  cnAt. 
Notice  that  the  updated  tracked  points  are  no  longer  the  true  closest  points.  Finally,  the 
Extension  Step  resets  each  nodal  value  to  be  a  true  closest  point. 


Implementation 

In  practice,  the  Initialization  Step  of  the  Closest  Point  Method  can  often  be  handled  ana¬ 
lytically  in  simple  problems.  More  complicated  wavefronts  can  be  treated  using  fast  tree- 
based  algorithms.  Implementation  of  the  evolution  Step  is  also  straightforward  because 
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each  tracked  point  is  just  updated  according  to  TP(x)new  =  XP{x)OTl9mal  +  cn&t.  The 
final  Extension  Step  is  more  complicated  and  is  typically  divided  into  two  parts,  a  search 
step  and  an  interpolation  step. 

In  the  search  step,  the  updated  value  fro  a  node  is  taken  to  be  the  closest  of  all  tracked 
points  (localization  of  this  step  is  possible).  This  gives  an  improved  approximation  of  the 
closest  point  representation.  Unfortunately,  this  process  cannot  create  any  new  tracked 
points  so  diverging  wavefronts  will  lose  resolution.  Thus,  a  second  interpolation  step  is 
needed  in  order  to  maintain  a  uniform  representation. 

Steinhoff,  Fan  and  Wang  carry  out  this  interpolation  by  averaging  over  nearby  nodes. 
This  very  simple  approach  is  effective  for  a  variety  of  interesting  problems,  but  it  can 
produce  spurious  wavefronts  in  certain  cases  and  is  low  order  accurate.  For  these  reasons, 
we  consider  a  higher  order  interpolation  based  on  nearby  neighbors.  These  neighbors  (call 
them  y  and  z)  are  chosen  so  that  x,y  and  z  are  collinear  and  roughly  parallel  to  the 
interface.  If  the  tracked  points  for  x  and  y  are  distinct  and  lie  on  the  same  smooth  curve, 
then  an  improved  estimate  for  the  closest  point  to  x  can  be  generated  using  the  nodal 
values  at  x  and  y.  Similarly,  an  improved  closest  point  estimate  can  be  attempted  using 
the  nodal  values  at  a;  and  z.  The  closest  of  these  two  results  to  x  is  taken  to  be  the  updated 
nodal  value. 

Notice  that  this  Extension  Step  does  not  yield  a  true  closest  point  representation.  How¬ 
ever,  closest  point  values  are  expected  near  the  interface.  Furthermore,  this  extension  has 
the  useful  property  that  every  nodal  value  represents  some  tracked  point  on  the  interface. 


First  Arrival  Times 

Unfortunately,  the  Closest  Point  Method  can  produce  gaps  in  the  interface.  We  now  discuss 
a  method  which  gives  a  much  more  uniform  representation  of  the  interface. 

The  formation  of  gaps  for  the  Closest  Point  Method  is  most  easily  understood  by  con¬ 
sidering  how  swallowtails  are  represented.  When  the  swallowtail  is  small,  nodal  values 
over-represent  corners.  Since  large  regions  are  used  to  represent  corner  points,  few  grid 
points  are  available  to  represent  the  end  of  the  swallowtail.  This  uneven  treatment  leads 
to  gaps  in  the  interface  which  propagate  and  grow. 

To  obtain  an  improved  result,  a  more  uniform  representation  is  needed.  For  example, 
we  can  set  each  nodal  value  to  be  the  point  on  the  interface  with  the  minimal  arrival 
time  rather  than  the  minimal  distance.  In  a  homogeneous  Medium,  without  reflection, 
this  just  means  that  each  nodal  value  is  set  equal  to  the  closest  point  on  the  interface 
that  propagates  directly  to  or  from  the  node.  Using  this  representation,  we  find  that 
redundancy  is  largely  eliminated,  and  a  greatly  improved  approximation  of  the  swallowtail 
is  obtained.  Unfortunately,  arrival  times  are  often  difficult  and  expensive  to  evaluate  in 
the  variable  index  of  refraction  case  or  when  reflections  occur.  Furthermore,  even  in  a 
homogeneous  medium,  this  approach  requires  a  more  intricate  search  step  nodal  values  can 
only  be  updated  when  a  nearby  tracked  point  travels  directly  towards  the  node  (which 
rarely  occurs). 

Of  course,  we  are  not  limited  to  representations  that  minimize  distance  or  arrival  times 
-  a  minimization  based  on  some  combination  of  distance  and  direction  of  motion  can  also 
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be  carried  out.  A  particularly  interesting  choice  arises  when  nodal  values  are  set  equal  to 
the  ”  Minimizing  Point” 

MP(x)  =  min  7|(x  -  y)  •  nx(y)|  +  ||x  -  y||2  (1) 

yelnterface 

for  7  >  0  since  then  a  good  agreement  with  the  minimal  arrival  time  representation  is 
found  near  the  interface.  Using  this  idea  leads  to  the  following  modification  of  the  Closest 
Point  Method: 


The  Arrival  Time  Method 

Initialize.  For  each  point  x:  Set  the  tracked  point  TP(x)  equal  to  the  minimizing  point 
MP(x)  on  the  initial  surface  r0  for  the  minimization  in  Eq.  (1).  Set  n  equal  to  the  normal 
at  MP(x). 

Repeat  for  all  steps: 

1.  Evolve  the  tracked  point  TP(x)  according  to  the  local  dynamics  for  a  time  At  : 
TP(x)t  =  ch. 

2.  Extend  the  surface  representation  by  replacing  each  tracked  point  TP(x)  by  the  point 
MP(x)  on  the  updated  surface  T  =  {TP(x)\x  €  Rn}  that  minimizes  Eq.  (1).  Re¬ 
places  each  n  by  the  normal  at  the  updated  MP(x). 

End. 

This  simple  approach  gives  a  much  more  uniform  representation  and  naturally  treats 
the  prototype  swallowtail  problem. 

We  have  localized  the  method  and  applied  it  to  problems  involving  refraction  and  re¬ 
flection.  We  have  also  calculated  the  intensity  by  using  wave  front  curvature.  The  results 
were  very  good. 


The  Segment  Projection  Approach 

Another  approach  to  this  problem  was  taken  by  A.-K.  Tornberg,  O.  Runborg,and  B.  En- 
gquist.  This  is  based  on  the  “segment  projection  method”  where  (a)  each  interface  7  is 
given  as  a  union  of  curve  segments  7 j.  (b)  The  segments  are  chosen  so  they  can  be  rep¬ 
resented  by  a  function  of  one  coordinate  variable,  (c)  The  domain  of  the  independent 
variables  of  these  functions  are  projections  of  the  segments  onto  the  coordinate  axes. 

In  addition  to  discretization  of  these  functions  one  needs  information  about  the  connec¬ 
tivity  of  the  segments.  For  each  part  of  an  x-segment,  one  needs  to  define  which  y-segment 
it  corresponds  to,  and  vice  versa. 

Each  of  these  functions  is  updated  according  to  a  simple  grid  based  advection  equation. 

This  technique  has  been  very  successfully  applied  to  a  host  of  variable  index  of  refraction, 
two  dimensional  problems.  Segments  cross  each  other  without  merging.  Swallowtails  and 
other  singularities  can  be  handled  by  the  use  of  phase  plane  variables. 
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